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ABSTRACT 


Using the Keck Planet Imager and Characterizer (KPIC), we obtained high-resolution (R~35,000) 
K-band spectra of the four planets orbiting HR 8799. We clearly detected H20 and CO in the 
atmospheres of HR 8799 c, d, and e, and tentatively detected a combination of CO and H20 in b. 
'These are the most challenging directly imaged exoplanets that have been observed at high spectral 
resolution to date when considering both their angular separations and flux ratios. We developed a 
forward modeling framework that allows us to jointly fit the spectra of the planets and the diffracted 
starlight simultaneously in a likelihood-based approach and obtained posterior probabilities on their 
effective temperatures, surface gravities, radial velocities, and spins. We measured vsin(i) values of 


Corresponding author: Jason J. Wang 
jwang4@caltech.edu 


2 WANG ET AL. 


10.1133 km/s for HR 8799 d and 15.0*27 km/s for HR 8799 e, and placed an upper limit of < 14 km/s 
of HR 8799 c. Under two different assumptions of their obliquities, we found tentative evidence that 
rotation velocity is anti-correlated with companion mass, which could indicate that magnetic braking 
with a circumplanetary disk at early times is less efficient at spinning down lower mass planets. 


Keywords: Exoplanet atmospheres (487), Exoplanet formation (492), High angular resolution (2167), 


High resolution spectroscopy (2096) 


1. INTRODUCTION 


In the past two decades, direct imaging has discov- 
ered several dozen substellar companions with masses 
from 1-70 Mjyp with orbital separations from ~3 au 
out to ~1000 au. (for a review, see Bowler 2016). The 
occurrence rates of giant planets and brown dwarfs be- 
yond = 10 au have begun to show that multiple forma- 
tion channels are responsible for the current population 
of imaged substellar companions (Nielsen et al. 2019; 
Vigan et al. 2020). Between 10-100 au, Nielsen et al. 
(2019) showed that giant planets between 5-13 Mjup 
have a higher occurrence rate compared to their brown 
dwarf counterparts (13-80 Mj,p) and preferentially oc- 
cur around higher mass stars, indicating the known exo- 
planet companions likely formed as the high-mass tail of 
planet formation through core accretion (Pollack et al. 
1996), whereas brown dwarf companions formed like bi- 
nary stars through gravitational instability (Boss 1998, 
2001). Vigan et al. (2020) found a similar dichotomy 
looking at the mass ratios between the companions and 
their host stars. T'hey inferred that companions around 
lower mass stars with mass ratios closer to unity formed 
like binary stars whereas more extreme mass ratio com- 
panions around more massive stars formed like planets. 

Since the directly imaged companions are amenable 
to spectroscopy and astrometric monitoring, we can 
use population-level characteristics beyond detection to 
study this population and understand how they formed. 
Bowler et al. (2020) reinforced the finding of multiple 
formation channels by showing evidence that giant plan- 
ets at wide separations (5-100 au) had an eccentricity 
distribution similar to that of close-in (« 1 au) giant 
planets, whereas the brown dwarf eccentricity distribu- 
tion resembled the stellar-binary population. Measure- 
ments of the spin of planetary mass companions have 
pointed to magnetic braking quickly slowing down the 
spin rate of planets after formation (Bryan et al. 2020a). 
While there has not yet been a population-level study of 
atmospheric compositions, compositional studies of in- 
dividual objects are able to contribute evidence to dis- 


* 51 Pegasi b Fellow 


cerning their formation channels (Konopacky et al. 2013; 
Barman et al. 2015; Gravity Collaboration et al. 2020; 
Molliére et al. 2020; Wilcomb et al. 2020). 

High-resolution spectroscopy of directly imaged com- 
panions allows us to characterize their orbits, spin, and 
compositions. The Doppler shift of molecular absorp- 
tion lines in the planetary atmosphere allows us to mea- 
sure the radial velocity of the planet, which is useful to 
break degeneracies between orbital inclination and ec- 
centricity for companions with limited orbital coverage 
(Snellen et al. 2014; Schwarz et al. 2016; Ruffio et al. 
2019). The rotational broadening of these absorption 
lines allows for a direct measurement of planetary spin 
(Snellen et al. 2014; Schwarz et al. 2016; Bryan et al. 
2018). The detection of molecular signatures through 
cross-correlation methods takes advantage of the fact the 
planet and star have different spectral features, enables 
the detection of trace molecular species, and allows for 
the inference of planetary composition (Konopacky et al. 
2013; Barman et al. 2015; Brogi & Line 2019; Wilcomb 
et al. 2020). 

However, up until now, slit spectrographs assisted 
with adaptive-optics (AO) have been observationally 
limited to bright companions that are at relatively large 
angular separations from their host stars. For planetary- 
mass companions within one arcsecond, only 8 Pic b 
and HR 8799 c have been detected at a spectral res- 
olution R > 10,000 (Snellen et al. 2014; Wang et al. 
2018a). The main difficulty is that the bright glare of 
the host star overwhelms the signal of the planet. The 
glare of the star can be mitigated by combining high- 
contrast imaging techniques with high-resolution spec- 
troscopy (Snellen et al. 2015). The combinations of these 
two techniques, which we term high dispersion corona- 
graph (HDC), drives the design of the Keck Planet Im- 
ager and Characterizer (KPIC; Mawet et al. 2017a; Jo- 
vanovic et al. 2019, Delorme et al. submitted). KPIC 
combines starlight suppression using the Keck AO sys- 
tem and single mode fibers with the NIRSPEC high- 
resolution spectrograph to enable high-resolution spec- 
troscopy of fainter and closer-in directly imaged plan- 
ets. Similar instrument designs are also being pursued 
by Subaru/REACH (Jovanovic et al. 2017; Kotani et al. 
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2020) and VLT/HiRISE (Vigan et al. 2018; Otten et al. 
2021). 

In this paper, we present the first science demonstra- 
tion of KPIC with observations of the four planets orbit- 
ing HR 8799. HR 8799 is a notable planetary system, as 
it is the only known system with four directly imaged ex- 
oplanets (Marois et al. 2008, 2010). The four planets are 
either near or in mean-motion resonance, and dynamical 
modeling of their orbits have constrained their masses 
to be 7.2 + 0.7 Mjup for the inner three planets and 
5.8 0.5 MJup for planet b (Wang et al. 2018b, also see 
Gozdziewski & Migaszewski 2020). Since their discov- 
ery, the planets have been characterized extensively in 
the 1-5 um range using broadband photometry and low- 
and medium-resolution spectroscopy (Bowler et al. 2010; 
Barman et al. 2011; Konopacky et al. 2013; Currie et al. 
2014; Ingraham et al. 2014; Skemer et al. 2014; Zurlo 
et al. 2016; Bonnefoy et al. 2016; Greenbaum et al. 2018; 
Molliére et al. 2020). At medium-resolution, the indi- 
vidual molecular lines begin to be resolved spectrally, 
allowing for detection of molecular signatures as well as 
the radial velocity of the planet (Barman et al. 2011; 
Konopacky et al. 2013; Ruffio et al. 2019). Measure- 
ments of elemental abundances in these planetary at- 
mospheres have shown deviations from the stellar abun- 
dances, which have been interpreted to mean that these 
planets formed via core accretion rather than gravita- 
tional instability (Konopacky et al. 2013; Barman et al. 
2015; Lavie et al. 2017; Molliére et al. 2020). 

Section 2 details the observations of the HR 8799 plan- 
ets made with KPIC. Section 3 describes the initial data 
reduction steps. The detection of molecular features 
through cross correlation as well as fitting atmospheric 
models of rotating planets directly to the data is dis- 
cussed in Section 4. We obtained the first spin measure- 
ments for these exoplanets, and we put them, as well 
as our measured orbital velocities and bulk atmospheric 
properties, in context in Section 5. We summarize our 
work and discuss future avenues both in obtaining better 
data and utilizing better models to study these planets 
in Section 6. 


2. OBSERVATIONS 
2.1. Instrument Description 


The Keck Planet Imager and Characterizer (KPIC) 
consists of a series of upgrades for the Keck II AO 
system and its two facility instruments: the NIRC2 
infrared imager and the upgraded NIRSPEC infrared 
high-resolution spectrograph (McLean et al. 1998; Mar- 
tin et al. 2018; López et al. 2020). As part of this project, 
instrument upgrades included a new vortex coronagraph 
for NIRC2 (Vargas Catalán et al. 2016; Serabyn et al. 


2017) and an infrared pyramid wavefront sensor operat- 
ing in H band for the AO system (Bond et al. 2020). Par- 
ticularly relevant for high resolution spectroscopy of di- 
rectly imaged exoplanets, a fiber injection unit (FIU) fol- 
lowing the concept presented in Mawet et al. (2017a) was 
deployed in 2018. The FIU benefits from the NIRSPEC 
detector upgrade that allows KPIC to reach R ~ 35,000 
(Martin et al. 2014, 2018; López et al. 2020). We point 
the reader to Delorme et al. (submitted) for detailed 
information about the instrumentation. 

Here, we provide a brief summary of the relevant in- 
strumentation. The pyramid wavefront sensor drives the 
facility deformable mirrors in the AO system to compen- 
sate for atmospheric turbulence. In addition to imaging 
exoplanets using NIRC2, KPIC can also send the K- 
band light of the system to the FIU to spectroscopically 
characterize these planets. Located after the AO sys- 
tem, the FIU steers the light of a planet into one of five 
single mode fibers represented by circles on Figure 1. 
These fibers are connected to NIRSPEC to spectrally 
disperse the light injected into the fibers (see right panel 
of Figure 1). Since the star is bright, diffracted starlight 
(stellar “speckles” ) leaks into all of the fibers, with the 
amount of starlight in each fiber depending on the an- 
gular distance between the fiber position on the sky and 
the star. In Mawet et al. (2017b), we demonstrated 
that the use of single-mode fibers provide the follow- 
ing advantages: a well-defined Gaussian-like line spread 
function (LSF) which is independent in shape to input 
wavefront aberrations, and, on average, 3x suppression 
of the underlying stellar speckle flux at the location of 
the fiber. Note that the current KPIC FIU does not 
utilize a coronagraph to suppress diffracted starlight, so 
the utilization of adaptive optics with the single mode 
fibers is the main starlight suppression mechanism. 

The fiber extraction unit (FEU) reimages the light 
from the single mode fibers onto the NIRSPEC spectro- 
graph slit, where the light is then dispersed inside the 
spectrograph. The resulting NIRSPEC data has up to 
five fibers illuminated with light in each echelle order 
(only four fibers were imaged onto the detector for the 
HR 8799 observations). Each fiber roughly subtends an 
angular diameter of ~ 50 mas. The projected separation 
between two consecutive fibers is ~800 mas for the four 
fibers shown. Because the slit sees the FEU and not 
the sky, the slit background is dominated primarily by 
the thermal emission from the warm room-temperature 
optics of the FEU, with a lesser contribution from the 
thermal instrument background of NIRSPEC. The sky 
background only appears in the fibers, but is generally 
well below the thermal background of the FEU, with 
only weak OH line emission seen in long exposures. 
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Fiber Positions 


NIRSPEC 


Figure 1. Schematic of KPIC observations of the HR 8799 planets. On the left are the locations of the fibers relative to a 
Keck/NIRC2 PSF of the star for an observing sequence on HR 8799 c. The red circle marks the location of the fiber centered 
on the planet, while the yellow circles mark the positions of the three other fibers that were in use. A fifth fiber is not shown 
in this image as it was blocked by the NIRSPEC slit and was not used. Note that the NIRC2 image was not taken during the 
observing sequence and is just to serve as an example to show the dynamic range of the stellar halo. On the right is a small 
portion (~360x370 pixels) of the NIRSPEC detector image in this configuration. A portion of two echelle orders are shown. In 
each order, four spectral traces are seen corresponding to the placement of the four fibers shown on the left image. 


2.2. HR 8799 Observations 


We listed the epochs of observations for each planet 
in Table 1. Each planet was observed with a similar 
observing sequence. NIRSPEC was set up to use the 
“Thin” filter, a thin piece of clear PK50 glass which 
blocks light at wavelengths longer than ~ 2.5 um, along 
with a custom pupil stop designed for the KPIC FEU 
and the NIRSPEC 0.0679”x1.13” slit to maximize the 
signal to noise of the light from each fiber relative to 
instrument background. The NIRSPEC echelle grating 
was set to 63.0? and the cross disperser was set to 35.76? 
to allow us to obtain nine spectral orders ranging from 
approximately 1.94 to 2.49 um. 

We started each observing sequence by placing the 
host star on each fiber and took at least one exposure 
(30-60 s) of the host star in all four fibers for data cali- 
bration purposes. We then designated one of the fibers 
as the primary science fiber, based on the fiber that had 
the best throughput in daytime testing. We used the 
tip/tilt mirror in the FIU to offset the star from the fiber 
bundle such that the planet of interest is placed on the 
primary science fiber. The science fiber received a com- 
bination of planet light plus diffracted residual starlight. 
The offset amplitude and direction was computed us- 
ing the whereistheplanet! orbit prediction tool (Wang 
et al. 2021) based on the dynamically stable and copla- 
nar orbit fit from Wang et al. (2018b). From prelimi- 
nary instrument characterization efforts, we estimated 
that the offset accuracy is 10 mas, which corresponds 
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to « 1096 loss in throughput due to fiber misalignment. 
Although the fibers in the bundle were fixed in a linear 
arrangement relative to one another, the adaptive optics 
field rotator (K-mirror) was used to rotate the field-of- 
view such that the star was coupled into as many of the 
other fibers as possible so that we could obtain simulta- 
neous stellar spectra. 

We then took 600 s exposures with NIRSPEC using 
the MCDS-16 detector readout mode in this configu- 
ration where one fiber had the light of the planet and 
at least two other fibers were transporting significant 
amounts of starlight.This integration time was chosen 
to be long enough so that read noise is negligible (ther- 
mal background noise ~10x larger), but short enough 
to tolerate any bad frames due to technical issues such 
as the AO loops opening. Every hour or so, we placed 
the host star on each science fiber and took a short ex- 
posure (30-60 s) for calibration. The open-shutter time 
obtained on each planet is listed in Table 1. 

Using the short exposures on the star, we calculated 
the end-to-end throughput from the top of the atmo- 
sphere to the detector. We reported the end-to-end 
throughput measured between 2.29-2.34 um in Table 1 
as a metric that combines both instrument performance 
and atmospheric conditions during the observations. We 
note that this wavelength range is not the best per- 
forming wavelength, as Delorme et al. (submitted) 
showed that a throughput of over 396 has been achieved 
at shorter wavelengths. However, our data analysis in 
the following sections is focused around this wavelength 
range since it coincides with the CO bandhead in K- 
band, so this is the most relevant throughput metric. 
Overall, we see that the conditions between the four 
nights are pretty comparable, with 2020-07-01 having 


HR 8799 PLANETS WITH KPIC 5 


Table 1. HR 8799 Planet Observations 


Date Target Integration Time (min) | Airmass | Throughput (%) 
2020-07-01 | HR 8799 c 230 1.0 - 1.8 1.7 - 1.9 
2020-07-02 | HR 8799 d 230 1.0- 1.7 1.9-2.4 
2020-07-03 | HR 8799 e 130 1.1- 2.0 2.2- 2.4 
2020-09-29 | HR 8799 b 160 1.0-1.1 1.8- 2.4 


slightly poorer performance due to issues in fiber injec- 
tion that were fixed after the night. 


2.3. Calibration Data 


In addition to the spectra of the star we obtained dur- 
ing the observing sequence, we also observed an M giant 
(HIP 81497) for wavelength solution calibration and a 
telluric standard star for the wavelength calibrator (ups 
Her) each night. We took five 1.5 s exposures of HIP 
81497 and three 30 s exposures of ups Her on-axis in 
each of the four fibers. After the night, we took thermal 
background frames looking at the FEU with no light in- 
jected into the fibers at each of the exposure times used 
to model the thermal background of the instrument. 


3. DATA REDUCTION 
3.1. Raw Data Reduction 


The process of going from the original detector images 
to extracted 1D spectra was the same for all data re- 
gardless of the object being observed. First, the images 
were background subtracted using the instrument ther- 
mal background frames taken when no light was being 
injected into the fibers during the daytime. The FEU 
typically was at a different temperature during the day 
so these thermal frames did not perfectly subtract the 
data and leave some residual background which we mod- 
eled during the extraction step. Thermal images taken 
during the observing sequence would have provided bet- 
ter background subtraction (i.e., nodding), but were not 
acquired as we had not developed an efficient way to nod 
the planet light on the detector. 

For each night of observation, the trace of each of the 
four science fibers in each of the nine orders was deter- 
mined by using the data on the telluric standard star ups 
Her. As the point spread function (PSF) of monochro- 
matic light coming from a single mode fiber is nearly 
a 2D Gaussian, we just needed to measure the position 
and standard deviation of the PSF at each wavelength. 
'To do that, we fit a 1D Gaussian to the cross section of 
the trace in each column of each order to determine the 
position and standard deviation of the Gaussian PSF at 
that column. For each of the nine orders, we recorded 


the position and standard deviation of the Gaussian PSF 
for each of the four fibers in each column of the detector 
(2048 in total). To mitigate measurement noise and bi- 
ases from telluric lines, we smoothed the measurements 
by fitting a cubic spline to the measured positions and 
standard deviations. These PSF standard deviations 
will also be used for estimating the line spread function 
(LSF) width in sections 4.2 and 4.3. 

For each exposure, we then extracted the 1D spectra 
of each of the four fibers in each column of each order. 
Due to the imperfect background subtraction, we used 
pixels that are at least 5 pixels away from the center of 
any fiber to estimate the residual background level in 
each column and subtracted the median of these pixels 
from every pixel in the column. Then, for each fiber, 
we used optimal extraction to measure the flux using 
a 1D Gaussian profile as the PSF with the positions 
and standard deviations fixed to the values measured 
on the telluric standard star (Horne 1986). The total 
integrated flux of the 1D Gaussian is then the flux we 
extracted for that fiber in that column. The uncertainty 
in the optimal extraction was used as the uncertainty in 
the flux measurement (Horne 1986). 


3.2. Wavelength Calibration 


We observed the M-giant HIP 81497 in each fiber 
each night to determine the wavelength solution from 
the stellar and telluric spectral lines. The wavelength 
solution was modeled as a spline using 6 nodes per or- 
der (i.e., piecewise-3'4 order polynomials). We modeled 
the data as the multiplication of a stellar spectrum, a 
telluric transmission term, and the transmission of the 
telescope and the instrument, complemented with an ad- 
ditive linear background term. The star HIP 81497, with 
a M2.5III spectral type, was chosen from the Gaia RV 
standard catalog (Soubiran et al. 2018). It was mod- 
eled with a PHOENIX stellar spectrum (Husser et al. 
2013) assuming a temperature of 3600 K, surface grav- 
ity of log(g) — 1, solar metallicity, and a fixed known 
RV of —55.567 + 0.0011 km/s. The telluric transmission 
of the Earth's atmosphere was modeled from a linearly 
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interpolated grid of 25 ATRAN? models (Lord 1992) 
over water vapor overburden (500, 1000, 5000, 10000, 
and 20000 um) and zenith angle (0, 25, 45, 68, and 89 
degrees). The water vapor overburden and the zenith 
angle were fit as nuisance parameters. The transmission 
of the telescope and instrument varies with wavelength 
primarily due to the efficiency of the spectrograph to dis- 
perse light at each wavelength (i.e., blaze function). To 
model the spectrally dependent transmission, we used 
a piecewise-linear function that divided each order in 5 
pieces. The best fit wavelength solution was found us- 
ing the Nelder-Mead optimization implemented in the 
scipy.optimize.minimize routine that jointly fit for 
the wavelength solution, telluric parameters, and the 
instrument and telescope transmission (Virtanen et al. 
2020). The search was initialized around the minimum 
of a coarse grid search to avoid local minima. This tech- 
nique was found to be precise to the 0.1 km/s level (Mor- 
ris et al. 2020), which will be sufficient for the following 
data reduction. Morris et al. (2020) also found that 
the wavelength solution was stable within a night to the 
same level of precision as long as optics inside the spec- 
trograph were not moved. 


4. FITTING HIGH-RESOLUTION SPECTRA 
4.1. Forward Modeling the Planet Spectra 


In this work, we build up a full forward model of the 
spectrum we recorded, which consists of both planet plus 
star light and is similar to the approach from Ruffio 
et al. (2019) for fitting medium-resolution integral field 
spectroscopy data of imaged exoplanets. The signal ob- 
tained from the fiber placed on one of the planets can 
be deconstructed into the following components: 


Dy(A) = a; (4)T (A) Prsr(A)+as(A)T (A) Sis p (A)9-n(A). 

(1) 
Here, D, is the signal from the planet fiber, T' is the 
transmission of the optical system (atmosphere, tele- 
scope, and instrument) excluding fiber coupling effi- 
ciency, Pr spr is the spectrum from the planet after it has 
been convolved by the instrumental LSF, a, is a scal- 
ing term for the planet brightness accounting for fiber 
coupling efficiency, Srsp is the spectrum from the star 
after it is convolved by the instrumental LSF (this will 
be modeled by the empirical spectra of the star), o, is 
a scaling term for the brightness of stellar speckles ac- 
counting for its fiber coupling efficiency, and n is the 
noise (instrumental thermal background noise generally 
dominates). For simplicity, we will drop the (A) notation 


? https:/ /atran.arc.nasa.gov/cgi-bin/atran/atran.cgi 


in the rest of the paper, but we note that all of these pa- 
rameters will remain wavelength-dependent implicitly. 
'The transmission of the optical system was calculated 
using on-axis observations of the star HR. 8799. In these 
observations, the signal D, can simply be written as: 


D,-—TSpysrp--n (2) 


We approximated $rpsp using a model PHOENIX spec- 
trum with an effective temperature of 7200 K (Husser 
et al. 2013). Given that HR 8799 is a FO star (Gray 
et al. 2003), it has nearly no spectral lines in K-band, 
which mitigates any errors due to an imperfect stellar 
spectrum. Additionally, we did not fit the data near the 
Bry line in our analysis. The signal-to-noise ratio (SNR) 
per spectral channel of the stellar spectra is 400, so n 
was assumed to be negligible. The transmission is then 
simply obtained by dividing the data by the model of 
the star: T = D,/Srsp. We used T solely to model the 
transmission of the planet light from the top of the atmo- 
spheres to the detector. To model the stellar light that 
leaked into the planet fiber, we simply used the on-axis 
observations of the star, Ds, to account for T(A) x S(A) 
simultaneously. This has the added benefit that the ex- 
act shape of the LSF only matters for the planet signal, 
and is not used in fitting the stellar spectrum. 

'Thus, we constructed a model for the data obtained 
on the planet fiber (Mp) as: 


My = pT Prsr + asDs (3) 


where we can use planetary atmosphere models for Pr sF 
to find the model atmosphere parameters that best fit 
the data. 

The terms a, and as vary slowly as a function of wave- 
length. The wavelength dependence of a, is dominated 
by differential atmospheric refraction changing the ap- 
parent sky position of the planet, and this effect changes 
slowly as a function of wavelength. The use of an at- 
mospheric dispersion corrector in KPIC Phase II will 
mitigate this effect on a, (Wang et al. 2020a; Jovanovic 
et al. 2020). We did not find chromatic optical aberra- 
tions in the system to be measurable for inclusion in ap. 
On similar observations of a speckle field with a single 
mode fiber, Gravity Collaboration et al. (2020) showed 
that o, can be approximated across the entire K-band 
as a low-order polynomial. In some preliminary analysis 
to characterize a, and a, from observations of standard 
stars, we confirmed this and found that they change on 
~0.01 um scales (~400 spectral channels). This means 
that high-pass filtering the data would be adequate for 
continuum subtraction and will remove chromatic mod- 
ulations in the continuum due to a, and as. 
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To do this high-pass filter, we median filtered the spec- 
trum with a 200-pixel (~0.004 um) box for computing 
the moving median and subtracted the median-filtered 
spectrum off from the original spectrum. We note that 
this continuum subtraction does not account for the fact 
that o; and a, also induce chromatic modulations in the 
line depths, which we will ignore in this work as it is a 
smaller effect. In principle, this effect can be modeled 
in our KPIC data using a low-pass filter approach dis- 
cussed in Ruffio et al. (2019). 

We applied the high-pass filtering to both the data 
and model spectra. The model spectrum after high pass 
filtering can be written as: 


H[M,| ~ ?t|os T Prsr + &sDs}, (4) 


where H denotes high pass filtering, and &p and à, are 
wavelength independent versions of a, and as. They 
represent achromatic scaling terms for the mean plan- 
etary and stellar speckle flux levels in the data. The 
high-pass filter of the data would be written in the same 
way. 

In Section 4.2 (but not in Section 4.3), we made one 
further approximation assuming that the cross-talk be- 
tween the various components is negligible when high- 
pass filtering using a median-filter implementation, al- 
lowing us to break down the model to filter the individ- 
ual components: 


H[M,| y ay T Prsr| + as? [D,] (5) 


We pulled à; and à; outside of H since we have assumed 
they are constant values after we removed the slowly 
varying chromatic continuum with the high-pass filter. 
We note that a linear implementation of the high-pass 
filter such as the Fourier-based approach from Ruffio 
et al. (2019) would make Equation 5 exact. However, 
we found that a median-filter implementation modeled 
the continuum better than the Fourier-based ones we 
tried (simple cutoff in frequency or a Gaussian filter in 
frequency space). Removing the chromatic continuum 
of o, and o, was more important than ensuring perfect 
linearity in the high-pass filter. 


4.2. Detection of Molecules 


First, we assessed the detection of all four planets by 
looking for the signature of molecular absorption lines 
from their planetary atmospheres as is common to do in 
the literature (e.g., Konopacky et al. 2013; Ruffio et al. 
2019; Xuan et al. 2020). As shown in Ruffio (2019), cross 
correlation techniques estimate the maximum likelihood 
value for the planet flux as a function of radial velocity 
(RV) shift for a given planet template Prsp. As we have 


spatially resolved data that allows us to construct the 
planet spectrum and stellar spectrum separately, we per- 
formed a modified cross-correlation like in Ruffio (2019) 
where we estimated the maximum likelihood value for 
both the planet and star flux as a function of RV shift 
for a given planet template. For a given RV shift, we 
can rewrite Equation 5 in matrix form to solve for the 
planet and star flux to best match the high-pass-filtered 
data H[D>]: 


AD,] | = | HIT Prsr(RV)| H[Ds] (>) . (6) 
Here, the data on the left hand side is a column vec- 
tor with a length equal to the number of spectral chan- 
nels, Ny. The model on right hand side consists of a 
N, x 2 matrix where the first column is the model flux 
of the planet and the second column is the model flux of 
the star (which we modeled empirically with the on-axis 
stellar spectrum), and a column vector of length 2 corre- 
sponding to the flux of the planet and star. The two un- 
knowns are à, and à,, making it straightforward to cal- 
culate their maximum likelihood values through linear 
least-squares optimization techniques. The cross corre- 
lation function (CCF) for a given planet model Pr sF 
is then the calculated values of oj at every requested 
RV. Note that we recomputed Prsr at every RV shift 
considered. When combining data from multiple echelle 
orders, the data and model columns can simply be ex- 
tended so that NX is the total number of spectral chan- 
nels from all of the orders considered in the fit. While 
not included in this analysis, each data point can also 
be weighted by its respective error when estimating à, 
by using equation D6 of Ruffio et al. (2019) for a more 
optimal “matched filtering" detection metric. In the 
following analysis in this section, such error weighting 
made negligible changes to the CCF (« 1096), so we 
have presented values without using it for simplicity. 

For Prsrp, we generated molecular templates from the 
cloudless Sonora-Bobcat model set (Marley et al. 2018, 
Marley et al. submitted). Although this model lacked 
clouds, we found that the CCF signal from the cloudless 
Sonora-Bobcat CO+H20 template was within 1096 of 
the CCF signal from the cloudy BT-Settl models used 
in Section 4.3 that contained all molecular opacities, in- 
dicating that the strength of our molecular detections 
does not strongly depend on cloud assumptions. Note 
that molecular templates constructed from the BT-Settl 
models are not available. We used the same atmospheric 
model parameters for all four planets: an effective tem- 
perature of 1200 K and a surface gravity of log(g) — 4.0 
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Figure 2. Cross correlation functions with individual molecular species for all four HR 8799 planets using the technique 
described in Section 4.2. For each planet’s spectrum (one planet per row), the SNR of the CCF are plotted for molecules 
templates using CO, H20, CH4, H20-4-CO (one per column). The solid color line is the cross correlation function of each 
molecular template with the data on the planet. The solid gray lines are the cross correlation function of each template with 
extracted spectra that contained no planet signal. The vertical dashed line marks the approximate systemic radial velocity as 
calculated by Ruffio et al. (2019). HR 8799 c, d, e have strong detections of CO and H2O individually, whereas HR 8799 b is 
only weakly detected when combining H2O and CO molecular templates. 


based on the latest values for HR 8799 e from Molliere 
et al. (2020). The detections were relatively insensitive 
to the exact choice of effective temperature and surface 
gravity: using a 1400 K model, which is representative 
of the scatter in effective temperature between different 
models (Greenbaum et al. 2018), changed the CCF sig- 
nal by 596?. Using the temperature structure and molec- 
ular abundance profile from this Sonora-Bobcat model 
and following the procedure in Morley et al. (2015), we 
post-processed the atmosphere model profile to compute 


3 Note, however, that the change in the likelihood of a model also 
depends on the number of data points so a small change in CCF 
signal can significantly alter the likelihood of a model (Brogi & 


Line 2019). 


emergent spectra at R = 200, 000, repeating the process 
to only consider the opacity of the molecule or molecules 
of interest each time. The model opacities used were 
from Freedman et al. (2008, 2014), which utilized the 
CH, line lists from Yurchenko & Tennyson (2014) and 
H50 line lists from Barber et al. (2006). The templates 
were then convolved to instrumental resolution assum- 
ing the spatial size of the fiber trace measured in Section 
3.1 is equivalent to the width of the LSF. 

We computed CCFs using CO, H20, CH4, and 
CO+H20 molecular templates. Spectra of each planet 
from echelle orders 31 (2.44-2.49 um), 32 (2.36-2.41 um), 
and 33 (2.29-2.34 um) were used, as these three orders 
had the best calibration. Including additional orders of 
the data or additional opacity sources in the model tem- 
plates did not significantly increase the CCF signal. For 
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each planet, we also computed CCFs for “noise” spectra 
that were taken contemporaneously. These noise spec- 
tra were either extracted from other fibers with similar 
amounts of stellar flux leaking into them as the science 
fiber but were not pointed at the planet, or from regions 
of the detector that imaged the slit without any fibers 
and thus were dominated by thermal background noise. 
The CCFs were computed with velocity offsets between 
-500 and 500 km/s from the Solar System barycenter. 
The CCFs for each planet were normalized by dividing 
by the standard deviation of the CCFs from all the noise 
spectra, resulting in CCF signal-to-noise (CCF SNR) 
functions. The CCF SNR functions for each planet and 
for each molecular species are plotted in Figure 2. We 
note that the CCFs of the noise spectra may not be 
Gaussian-distributed, so there maybe biases when quan- 
tifying the false alarm probabilities of these molecular 
detections. 

We found strong detections of CO for HR 8799 c, d, 
and e with SNR between 7-11. For these three plan- 
ets, we also have a strong detection of H20 with SNR 
> 5. For HR 8799 b, we have a weak detection with 
SNR z 3 when combining the CO and H20 templates, 
but there is no significant detection of any individual 
molecule alone. This is likely due to the fact HR 8799 b 
is about twice as faint as the other planets and had one 
of the shortest exposure times. Longer exposure times 
are needed to obtain confident detections of molecular 
signatures in HR 8799 b. For all four planets we did 
not detect CH; at a significant level. The non-detection 
of CH, despite the strong CO and H20 detections are 
consistent with molecular detections of these planets at 
medium resolution with OSIRIS (Konopacky et al. 2013; 
Barman et al. 2015; Petit dit de la Roche et al. 2018, 
Ruffio et al. submitted). Our KPIC detections are the 
first detections of HR. 8799 b, d, and e at high spectral 
resolution (R > 10,000), where many molecular absorp- 
tion features start becoming spectrally resolved. Our 
6c detection of H20 in HR 8799 c in K band is bet- 
ter than the previous 4.60 L-band detection of H20 by 
NIRSPAO that used 3.5x more integration time (Wang 
et al. 2018a). 


4.3. Bayesian Inference of Planetary Parameters 


We put our forward modeling approach into a 
Bayesian inference framework to fit our extracted spec- 
tra directly, retrieve atmospheric parameters, and assess 
how complex of a model is required to fit the data. Our 
likelihood-based framework is similar to the one devel- 
oped by Ruffio et al. (2019) for medium resolution spec- 
troscopy of the HR 8799 planets, although we did not 
analytically marginalize over any of the parameters in 


our likelihood. Unlike Brogi & Line (2019), we did not 
use the cross correlation function in the likelihood, as we 
did not assume that each spectral channel has the same 
noise level. Our method is similar to the method used 
by Gibson et al. (2020) to characterize the Fe absorp- 
tion on an ultra-hot Jupiter, except we simultaneously 
fit the star and planet together, which minimizes over- 
fitting of the planet signal when subtracting off the star 
using techniques such as principal component analysis 
(Pueyo 2016; Ruffio et al. 2019). We further note that 
the relative flux ratio between the planet and the stel- 
lar components in the data are much less extreme in our 
case, which likely makes it easier to fit both components 
to the data simultaneously. 

We were interested in constraining the atmospheric 
and bulk parameters of the planets. For the planet 
spectrum, we fit the B T-Settl-CIFIST model grid to the 
data, varying both effective temperature (Ter) and sur- 
face gravity (log(g) in cgs units) (Allard et al. 2012). 
We chose the BT-Settl models as it was the only pub- 
licly available grid of models that are available at a spec- 
tral resolution R > 35,000 and includes clouds, which 
have been shown to be important shaping these planets' 
spectra (e.g., Molliére et al. 2020). In particular for our 
high resolution spectra, the depths of molecular absorp- 
tion lines change when clouds are included (Hood et al. 
2020). Rather than just stepping through RV shifts as 
in Section 4.2, we fit for the planetary radial velocities 
relative to the Solar System barycenter. Note that the 
stellar radial velocity is not well determined making rel- 
ative RV measurements challenging (Ruffio et al. 2019). 
Similar to Section 4.2, we can fit for the planet fluxes, 
which controls the depth of the planetary atmosphere 
lines compared to the stellar and telluric line depths. 
This can also be directly translated to K-band flux ra- 
tios of the planets which we can check against litera- 
ture photometric measurements. We also fit for the ro- 
tational broadening of the planetary spectra (vsin(i)) 
by using the fastRotBroad function in PyAstronomy to 
broaden our planet atmosphere models (Czesla et al. 
2019). 

We also improved on the stellar model by considering 
multiple D, exposures. In Section 4.2, we averaged all of 
the on-axis images of the star in time. Here, we fit a set 
of linear coefficients c; to optimally weigh the individual 
on-axis exposures of the star, D,;, to create the master 
stellar spectrum that best fits the data: 


D; = 2, eiDsil 2 Ci. (7) 


We note that this is similar to the LOCI technique in 
high-contrast imaging (Lafreniére et al. 2007), except 
that we optimized the coefficients while simultaneously 
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fitting the planet model to the data like what has been 
done in medium-resolution spectroscopy (Ruffio et al. 
2019). In our fits, we assumed the c; coefficients to be 
unchanged across orders, but the overall spectrum D, 
can be scaled by a different flux value to account for the 
chromaticity of the stellar speckle flux. 

Accurate estimation of errors is required for any ro- 
bust statistical analysis. Preliminary analysis of the 
data indicated that the residuals to the forward model 
fits are dominated by uncorrelated noise. However, we 
found that the amplitude of the uncorrelated noise was 
higher than what the formal extraction errors predicted. 
This may be due to an underestimation of extraction er- 
rors or due to unaccounted for noise terms by ~30%. Fu- 
ture work that performs more thorough analysis of the 
instrumental noise and data pipeline could help identify 
the source of this noise. For the purpose of this work, we 
simply fit for this excess uncorrelated noise, assuming it 
is Gaussian. Thus, the total error of, = Ofipe + Ohy 
where dpipe is the nominal extraction error from our 
pipeline and og; is the error term we fit for. We used a 
separate og, for each spectral order, but assumed that 
Ost is constant within an order. This seemed to be 
a suitable approximation based on our analysis of the 
residuals to the fit. 

In the analysis until now, we assumed that the width 
of the Gaussian trace of the fiber in the spatial dimen- 
sion could be used as the LSF in the dispersion di- 
mension. This assumption is not perfect as the NIR- 
SPEC spectrograph was designed with a difference in 
focal lengths in the spatial and dispersion directions by 
a factor of 1.13 (Robichaud et al. 1998). In prelimi- 
nary analysis of the OH sky lines which are unresolved 
at our resolution, we found that the LSF in the disper- 
sion direction is indeed 1.12 + 0.02 times wider than the 
width of the Gaussian profile measured in the spatial 
dimension. To conservatively account for any systemat- 
ics in this preliminary measurement, we will allow the 
LSF width to vary between 1.0-1.2 times the width we 
measured in the spatial direction in Section 3.1 (i.e., the 
aspect ratio of the 2D LSF). This corresponds to varying 
the resolution from ~35,000 to ~29,000. We note that 
since the stellar spectrum is built using empirical data, 
the LSF size only affects the broadening of the planetary 
model, Prsr. 

We defined the log-likelihood to be: 


mee) = -3 (AP + In no) 


2 OF ot 
(8) 
where we are summing over each spectral channel in the 
data. ^4[M,] was constructed using Equation 4. Note 
that implicit in M, are the parameters of the planet, the 


instrumental LSF, and the nuisance parameters of the 
stellar speckle spectrum. For each planet, we consid- 
ered three different models: a forward model that only 
contains the stellar speckle spectrum and no planet sig- 
nature (i.e., &p = 0 and all planet parameters are fixed) 
which we call “No Planet”, a model containing a planet 
with no discernable rotation (vsin(i) — 0 is fixed) which 
we call *No Rotation", and a rotating planet model 
where all planet parameters are allowed to vary which 
we call “Full Model”. 

For HR 8799 b, c, and d, we used echelle orders 31 
(2.44-2.49 um), 32 (2.36-2.41 um), and 33 (2.29-2.34 
um) to fit our model to the data just as was done in Sec- 
tion 4.2. We did not use orders 34-39 (1.95-2.29 um) as 
they either had poor wavelength calibration (uncertain- 
ties larger than a spectral channel), or had strong CO% 
telluric absorption that was difficult to forward model 
accurately. For HR. 8799 e, we only used orders 32 and 
33 as we found that the strong telluric absorption fea- 
tures in order 31 could not be fully modeled, with cor- 
related residuals of comparable amplitude to the planet 
signal. We omitted this order to mitigate the effect of 
residual telluric lines from biasing our fit of the planetary 
atmosphere. Future work to marginalize over localized 
telluric residuals (Czekala et al. 2015) or improvements 
in modeling the stellar speckle spectrum could make this 
order useful in the fit, but we chose to omit this order 
for now. 

For the free parameters in each model, we used the 
following priors. For the planet properties, an uniform 
prior on Tg between 1000 and 1800 K, an uniform prior 
on log(g) between 3.5 and 5.5, an uniform prior on the 
planetary radial velocity between -150 and 150 km/s, an 
uniform prior on v sin(i) in the Full Model between 0 and 
60 km/s, and an uniform prior on the planet flux from 0 
to 25 DN (data numbers). For each order, we included 
a parameter for the stellar speckle flux and nuisance pa- 
rameters to fit for systematics in the data: an uniform 
prior on the stellar speckles flux between 0 and 500 DN, 
an uniform prior on a residual background term between 
-10 and 10 DN due to imperfect background subtrac- 
tion, and a log-uniform prior on the og, between 0.1 
and 30 DN. In the Full Model and No Rotation models, 
we also included for each order a multiplicative factor 
to account for any broadening of the LSF beyond what 
we measured as the spatial width of the fiber trace with 
an uniform prior between 1.0 and 1.2. This term was 
not in the No Planet model because that model con- 
sisted solely of empirical data where we did not need 
to broaden anything to instrumental resolution. To fit 
the stellar speckle spectrum using a linear combination 
of on-axis stellar spectra, we used an uniform prior be- 
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tween 0 and 2 for each c; term. If we term N, as the 
number of stellar spectra, D, ;, used to compute D, and 
Norder as the number of spectral orders used in the fit, 
the No Planet model has 3Norder + Ns free parameters, 
the No Rotation model has 4+4Norder + N, parameters, 
and the Full Model has 5+4Norder + Ns free parameters. 
In the end, this resulted in ~20-30 free parameters, of 
which, we were mostly interested in the planetary pa- 
rameters. 

We sampled the posterior using the nested sam- 
pling algorithm (Skilling 2004, 2006) implemented in 
pymultinest, which allowed us to both perform param- 
eter estimation and compute the model evidence (Buch- 
ner et al. 2014). 'The evidence, P(D|M), is used to 
compute the Bayes factor, B, that can be used to assess 
the relative probability of model M4 compared to Mg. 
If P(Mi|D)/ P(M3|D) expresses the relative probability 
of the two models given the data, then, 

. P(D|Mi) | P(Mi|D)P(M3) 


B = Pw) PüspPauy — 


where P(M) is the prior probability of that model. As 
we assumed each model has equal prior probability, then 
B is equivalent to the relative probability of two mod- 
els. We used the Bayes factor to compare the simpler 
models to the Full Model to determine whether the ad- 
ditional free parameters are justified. Thus, our M» will 
always be the Full Model. We listed the estimates on the 
planet parameters and the Bayes factor of each model 
compared to the Full Model of that planet in Table 2. 
For the Full Model fits, we also plotted the posterior dis- 
tributions of the key planet parameters in Figure 3. For 
the Full Model fits, the strongest covariance is between 
Te and log(g) which we discuss in Section 5.1. We show 
corner plots of the posterior distributions from the Full 
Model fits in Appendix A. 


For each planet, we can “decisively” reject any model 
with a B that is more than 100 times smaller than the 
model with the highest B based on the interpretation of 
the Bayes factor suggested by Jeffreys (1983). In Table 
2, we see that the No Planet model is ruled out for c, d, 
and e, but remains 9% as likely as the Full Model for HR 
8799 b. This finding is consistent with the marginal 3o 
detection of b using template cross correlation in Section 
4.2, but we are now able to assign relative probabilities 
to the cases. The Bayes factor between a planet and no 
planet model offers an alternative way to determine de- 
tection significance rather than using CCF SNR, where 
the false alarm probability is unclear. Bayesian hypothe- 
sis testing methods have been used previously for source 


detection in cosmological datasets (Carvalho et al. 2009) 
and for exoplanet direct imaging (Golomb et al. 2019). 

The data is less definitive in distinguishing between 
No Rotation and Full Model. In all cases, the No Rota- 
tion has a > 1% probability compared to the Full Model. 
We know that it is unphysical to assume these planets 
have no spin, but the No Rotation model is a good ap- 
proximation for a planet with a spin that remains un- 
detectable. Thus, we interpret this result as the current 
data does not provide a definitive detection for rota- 
tional broadening. This could be due a low signal to 
noise detection like in the case of HR 8799 b, but it could 
also be the difficulty of measuring small vsin(i) values, 
especially given the near-face-on orbital configuration of 
the system. Still, HR 8799 d and e have a No Rotation 
model with B « 0.05, which Jeffreys (1983) interprets 
as “very strong" evidence for rotational broadening: the 
Full Planet model being 30 times more likely for d and 20 
times more likely of e. We note that vsin(;) values from 
the Full Model for both HR. 8799 d and e are inconsis- 
tent with 0 by > 30, but the B metric downweighs this 
to account for the addition of this free parameter that 
could cause overfitting. B provides a more straightfor- 
ward and more conservative assessment of the detection 
of rotational broadening rather than determining how 
far vsin(i) must deviate from 0 to be quantified as a de- 
tection, since the median vsin(?) will always be nonzero 
by construction. For HR 8799 b and c, we derive a 95% 
upper limit on the vsin(;) of 51 and 14 km/s, respec- 
tively, based on their marginalized 1D posteriors plotted 
in Figure 3. Thus, in this work, we were able to make 
the first measurements or constraints of the rotation ve- 
locities of the HR 8799 planets. In the future, higher 
SNR detections or detections at higher spectral resolu- 
tion would enable more definitive detections of rotation. 

To assess the quality of our fits, we plotted the best 
fitting set of parameters from the Full Model fit to one 
order of the data for all four planets in Figure 4. Visual 
inspection indicated our forward models fit the data ad- 
equately, with the residuals appearing to be dominated 
by uncorrelated noise. We verified this by computing 
the autocorrelation function (ACF) of the residuals and 
found that the ACF is well approximated by a delta 
function, with the wings of the ACF having an ampli- 
tude of <5% of the peak (see Appendix B). This con- 
firms that an uncorrelated noise model is sufficient, as 
we are likely dominated by thermal background noise of 
the instrument in our observations. 


5. DISCUSSION 


5.1. Atmospheric Parameters 
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Table 2. Model fits to HR 8799 KPIC data. For each parameter, the median value is listed, with the subscript and superscript 
values representing the range of the central 68% credible interval with equal probability above and below the median (the central 
95% credible interval is listed in parentheses). 


Planet Model Tor (K) log(g) RV (km/s) vsin(i) (km/s) Planet Flux (DN) B 
+278 .4(+362.3) +0.4(-40.7) +6.5(+12.5) T11.9(022.5) T1.9(13.6) 
b Full Model | 1423.3 5,26 3540) 48 ga 15  —19.6 710 149). 312 1 z( 224) 9.4 V 7(—3.0) 1.0 
: +275.4(4+339.1) 4-0.5(4-0.8) +3.2(+6.6) +1.4(+2.8) 
b No Rotation 1448.5 369.1(—399.4) 4T 7-11) —9.8 5 o(—16.3) — 3-4 33 2,0) 0.13 
b No Planet — — — — — 0.092 
126.2(182.2) +0.1(-F0.1) +0.8(+1.6) T3.8(07.3) +1.1(+4.0) 
c Full Model Dade uos Sa. aa ei gat as 1.0 
, --24.4(4-38.9) +0.1(-+0.1) +0.9(+1.5) +0.7(+1.5) 
c No Rotation 1474.4 36.3(—102.6) 5.4_0.2(—0.4) —12.4^ 1 ot—1.9) — 9.0 o 6513) 2.1 
c No Planet — — — — — 1.6 x 107 
+50.9(+105.6) X0.3(T0.4) +1. 1(42.2) +2.8(45.7) X23 44.8) 
d Full Model 1558.8 e: ane TA Mii TEE ee M8 ms a 1.0 
; +108.8(+150.5) +0.1(-+0.1) +1.4(+42.7) 1.8(+3.8 
d. No Rotation | 189. 1 Gea ey, leet = nut 0.032 
d No Planet — — — — — 3.1 x 1077? 
X57.0(T124.8) +0.3(-0.9) 31.2(1:2.5) X2.3(74.6) ¥4.5(47.7) 
e Full Model 145.6 55 ao ST cna -iagt aas) iugo prea 162 56 s 1.0 
; +124.3(+318.3) +0.7(+1.1) +1.1(+2.6) +1.7(+5.4 
e No Rotation 1323.07 37 G( 163.8) 4.3 06-08) —H.6 Vi 1) = 9.27 1 3(-2.3) 0.049 
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Figure 3. Posterior distributions of the main planetary parameters for each planet in the Full Model fit. In the bottom left 
panel, the gray shaded region represents the 6896 credible interval of the radial velocity of the host star relative to the Solar 
System barycenter. 
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Figure 4. The 1D extracted spectra (black) from the fiber placed on each planet and the best fitting forward model from echelle 
order 33. The forward model (blue) has been deconstructed into its two constituent parts: the stellar model (cyan) built from a 
linear combination of on-axis stellar spectra and the planet model generated from the BT-Settl atmospheric models (red). The 
residuals to the fit are plotted as gray circles, and appear to be dominated by uncorrelated noise. A zoomed-in version of this 
plot is available in Appendix B. 
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The detection of CO and H20 but not CH; in our 
high-resolution spectra is consistent with previous at- 
mospheric studies of the HR 8799 planets. Our CCF 
detections of HR 8799 c and d agree well with previous 
molecular cross-correlation analyses at lower resolution 
(Konopacky et al. 2013; Wang et al. 2018a, Ruffio et al. 
in prep.). Due to the weak detection of HR 8799 b in our 
limited observations, we were only able to marginally 
detect the combined signal of CO and H20 and did not 
have enough SNR to address previous disagreements on 
the amount of methane in its atmosphere (Barman et al. 
2015; Petit dit de la Roche et al. 2018). Longer inte- 
gration times and improvements to instrument perfor- 
mance will improve on these data. Regardless, this is the 
first time HR 8799 e is studied at high spectral resolu- 
tion, which demonstrates the high-contrast capabilities 
of KPIC. Previously, the highest resolution spectrum for 
HR 8799 e was the R ~ 500 GRAVITY spectrum (Grav- 
ity Collaboration et al. 2019). For HR 8799 e, we found 
detections of CO and H20 of similar strength as plan- 
ets c and e and a similar non-detection of CH4. This is 
consistent with the picture that the inner three planets 
have similar spectral signatures and luminosities based 
on lower resolution spectroscopy. 

We compared the planet fluxes measured in our for- 
ward model fits and reported in Table 2 to the ex- 
pected fluxes of the planets using the end-to-end sys- 
tem throughputs reported in Table 1, the exposure times 
of each frame, and the gain of the NIRSPEC detector 
(3.03 e-/DN; López et al. 2020). Using the photome- 
try of the planets in the SPHERE K2-band from Zurlo 
et al. (2016), we expected to measure a flux of 4, 13, 
16, and 17 DN for HR 8799 b, c, d, and e respectively. 
These values are mostly consistent with our measured 
fluxes, although it appears that our inferred fluxes for 
HR 8799 c, d, and e to be lower on average. This could 
be due to losses in flux due to a misaligned fiber when 
the blind offset was performed. It is not entirely clear 
this is the case, since we do not see a similar effect for 
HR 8799 b, which had the largest blind offset. Still, the 
9596 credible intervals of the fluxes are consistent with 
the literature photometry for all four planets. The fact 
the planet fluxes measured from our high resolution data 
agrees with broadband photometry increases our confi- 
dence that we have reliably separated the planet signal 
from the star. 

The bulk atmospheric properties (effective tempera- 
ture and surface gravity) we obtained from our forward 
model fits somewhat agree with previous work at lower 
resolutions. Our Teg between 1400 and 1700 K for HR 
8799 c and d agrees quite well with those obtained by 
Bonnefoy et al. (2016) and Greenbaum et al. (2018) in 


their BT-Settl fits, but HR 8799 e with a Tug between 
1200 and 1500 K is lower than those works by 300 K 
(only 0.0296 of our posterior have effective temperature 
> 1650 K). We note they obtained nonphysical radii 
in their BT-Settl fits which pinpoint to issues with fit- 
ting the BT-Settl grid to the HR 8799 spectral data (we 
did not perform absolute flux calibration so we did not 
constrain their radii). The log(g) values for HR 8799 
e agree well with Bonnefoy et al. (2016) and Green- 
baum et al. (2018). The posteriors for HR. 8799 c and 
d favor log(g) > 4 and peak at the upper bound of 
the prior (i.e., upper bound of model grid), which is at 
odds with those works which prefer log(g) « 4 (only 
0.0596 and 396 of our posteriors for HR. 8799 c and d, 
respectively, had log(g) « 4). When looking beyond the 
BT-Settl model fits, our effective temperatures are sys- 
tematically higher than the ~1100 K derived from other 
atmospheric models or predicted by evolutionary mod- 
els, and log(g) values for HR. 8799 c and d continue to 
remain higher than literature values between 3.5 and 4 
(Marois et al. 2008, 2010; Konopacky et al. 2013; Currie 
et al. 2014; Ingraham et al. 2014; Skemer et al. 2014; 
Bonnefoy et al. 2016; Greenbaum et al. 2018; Molliére 
et al. 2020). This indicates there may be some system- 
atic errors in the BT-Settl models, since both high and 
low resolution data both found higher 7;g values for the 
BT-Settl model than all other models. One reason could 
be that BT-Settl models use equilibrium chemistry, and 
there has been evidence for disequilibrium chemistry in 
the HR 8799 planets causing stronger CO lines than ex- 
pected and could make equilibrium chemistry models 
favor higher temperatures (Skemer et al. 2014; Molliére 
et al. 2020). However, we were not able to fit our data to 
other publicly available model grids with clouds as they 
are not generated at sufficiently high spectral resolution 
to match our data. 

HR 8799 c, d, and e are thought to have similar bulk 
atmospheric properties (e.g., luminosity, effective tem- 
perature, radius) based on broadband spectroscopy and 
photometry (e.g., Marois et al. 2010; Bonnefoy et al. 
2016). This appears to be at odds with our marginal- 
ized 1D posteriors, which prefer different Teg and log(g) 
values for HR 8799 e. Partially, this is due to a strong 
correlation between log(g) and Tig: lower Teg values in 
our fits also prefer a lower log(g) (left panel of Figure 
5). The contours for all three planets lie on nearly the 
same plane, suggesting that the differences in both Teg 
and log(g) may stem from additional parameters that 
we were not able to adjust in the model grid fits (e.g., 
cloud properties or disequilibrium chemistry). The need 
to adjust more parameters has been seen in B T-Settl fits 
to substellar companions where additional extinction pa- 
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Figure 5. Atmospheric constraints on the HR 8799 planets. On the left are 1D and 2D marginalized posteriors of Teg and 
log(g) for HR 8799 c (red), d (green), and e (yellow). The 2D plot shows the 68% (solid) 9596 (dot-dashed) credible region for 
each planet. A strong positive correlation between Tig and log(g) are observed. Black points and dotted lines correspond to 


representative models, which are plotted on the right. 


rameters have been used to augment the default cloud 
prescription (Marocco et al. 2014; Bonnefoy et al. 2016; 
Delorme et al. 2017; Ward-Duong et al. 2021). 

Looking more closely at why these fits to just our high- 
resolution spectra prefer these values, we plot how the 
BT-Settl spectra change as we change surface gravity 
and effective temperature in this correlated way in the 
right panel of Figure 5. We observed that the spectra are 
quite similar visually. The Teg = 1300 K, log(g) = 3.5 
model preferred by HR 8799 e has similar CO line depths 
in many of the lines as the Teg = 1650 K, log(g) = 5.5 
model preferred by HR 8799 d, whereas the interme- 
diate temperature models have the deepest CO lines. 
Since the BT-Settl models assumed solar abundances in 
carbon and oxygen, such changes in the H2O and CO 
line depths could also be due to differences in abun- 
dance, given the broadband spectral data indicates the 
planets should have similar temperatures. Performing 
abundance measurements on this high resolution data 
in the future will offer independent and sensitive mea- 
surements of the C/O ratios of the planets. 

While we did not observe that planetary RV and 
vsin(i) are strongly correlated to Tyg and log(g) (see 
Appendix A for corner plots), any inconsistencies in 
these bulk atmospheric properties could bias our in- 
ferred planetary RV and vsin(i) values. This should 
partially be mitigated by the fits having broad priors in 
Ter and log(g), which allowed us to marginalize over a 
large range of possible models in deriving 1D posteriors 
for planetary RV and vsin(i). However, there could be 
remaining systematic errors in the model that are not 
accounted for. To examine if there could be additional 
biases, we fixed Tyg and log(g) to the literature best- 
fit values for the BT-Settl grid from Greenbaum et al. 
(2018) for HR. 8799 c, d, and e, and reran the fits. We 
found that planetary RV changed by < lo and vsin(i) 


changed by < 1.50 for all three cases. We also assessed 
the sensitivity to cloud assumptions by rerunning the fits 
using the cloudless Sonora-Bobcat models instead, and 
found that both planetary RV and vsin(i) changed by 
« lo. These changes are comparable to our statistical 
errors, so any systematic biases should not significantly 
alter our results on RV and spin. We note that the B of 
these fits were < 107? compared to the Full Model, so 
cloudless models are not good fits to the data, even if the 
planetary RV and vsin(i) appear to be consistent. Any 
smaller tweaks to the spectra (e.g., changing the C/O 
ratio) should have even smaller effects on the inferred 
planetary RV and v sin(i). 


5.2. Orbital Velocity 


The planetary radial velocities we measured are not 
precise enough to add significant constraints on their 
orbits, but are consistent in amplitude and sign with 
previous work (Wang et al. 2018b; Ruffio et al. 2019). 
Using the dynamically stable solutions from Wang et al. 
(2018b) and the measurement that the position angle 
of the ascending node is less than 180 degrees (Ruffio 
et al. 2019), we expect an RV offset relative to the star 
of 1.9+ 0.2 km/s for planet b, 0.0+ 0.2 km/s for planet 
c, —2.7 E 0.3 km/s for planet d, and —2.1+0.3 km/s for 
planet e. We adopted the convention that negative RVs 
indicate the planet is moving towards us. The stellar 
radial velocity is somewhat uncertain with the latest es- 
timate being Eu km/s with respect to the Solar 
System barycenter (Ruffio et al. 2019). 

Combining these two terms together, we have 
predicted planetary RVs relative to the Solar Sys- 
tem  barycenter of —8.6*05 km/s for planet b, 
—10.5*05 km/s for planet c, -32t05 km/s for planet 
d, —12.6*02 km/s for planet e. Although within the 
95% credible interval of all our measurements, our val- 
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ues for HR 8799 b, c, and d are systematically more neg- 
ative by ~lo (~1 km/s). This could indicate an error 
in the absolute wavelength solution of our instrument, 
or systematics in computing the stellar radial velocity 
in previous analyses at the 1 km/s level. Overall, these 
results confirm the ability to measure planetary radial 
velocities with KPIC. 


5.3. Planetary Spin 


We used our measured vsin(i) values to compare the 
rotation velocities for HR. 8799 c, d, and e to the popu- 
lation of substellar companions and free-floating objects 
with measured spins. Although the spin measurements 
for HR. 8799 c and d are not definitive, the strong evi- 
dence for rotation still gave us confidence to assess what 
these measurements would imply about their rotational 
evolution. We excluded HR 8799 b in the analysis since 
its vsin(i) is unconstrained given the tentative detec- 
tion, and any preference in vsin(i) values may be spuri- 
ous (we still reported some numbers taking the posterior 
at face value). 

The bulk of the companion population was studied 
in Bryan et al. (20202). We also did not consider HD 
106906 b and AB Pic b in our subsequent analyses, since 
they only have tenuous rotation measurements (Zhou 
et al. 2019, 2020). We included the vsin(i) measured 
for GQ Lup b (Schwarz et al. 2016), the rotation pe- 
riod measured for HN Peg b (Zhou et al. 2016), and 
the rotation period of Jupiter (Dessler 1983) and Saturn 
(Helled et al. 2015) to extend the mass range of substel- 
lar companions considered. For the free-floating objects, 
we included the planetary mass objects from Bryan et al. 
(2020a) as well as variable brown dwarfs with photomet- 
rically measured rotation periods compiled in Vos et al. 
(2020) and supplemented with measurements from Tan- 
nock et al. (2021). For the variable brown dwarfs, we 
only used those with parallax measurements from the lit- 
erature (Faherty et al. 2012; Smart et al. 2013; Tinney 
et al. 2014; Liu et al. 2016; Theissen 2018; Gaia Col- 
laboration et al. 2018; Best et al. 2020; Gaia Collabora- 
tion et al. 2020) in order to convert their K-band fluxes 
(Cutri et al. 2003) to masses using evolutionary tracks 
(Baraffe et al. 2003). We note that this sample contains 
objects formed through a variety of formation mecha- 
nisms: free-floating objects likely are the low-mass tail 
of cloud fragmentation (e.g., Luhman 2012), while plan- 
ets like Jupiter and Saturn formed via core-accretion 
(Pollack et al. 1996). Thus, depending on the dominant 


processes that produced the spins we measured, these 
populations may or may not obey similar spin trends. 
However, Bryan et al. (2020a) recently argued for a sin- 
gle spin regulation mechanism for both free-floating ob- 
jects and companions. 

Converting v sin(i) measurements requires accounting 
for the inclination degeneracy, and for most of these ob- 
jects, including the HR 8799 planets, there are no em- 
pirical constraints on the spin axis inclination, and thus 
obliquity. While Vos et al. (2020) found a correlation 
between near-infrared color anomaly with spin axis in- 
clination for low-gravity brown dwarfs, this relation is 
not particularly constraining for the HR 8799 planets: 
they have (J — K) ~ 2.5 mag (Zurlo et al. 2016) that 
corresponds to a (J — K) color anomaly of ~ —0.1 mag 
(Liu et al. 2016), which lands in a region where incli- 
nations between 20? and 90? were measured (Vos et al. 
2020). Directly measuring the inclination of the spin 
axis typically requires a photometrically derived rota- 
tion period and a rotational broadening measurement, 
and an obliquity constraint needs the spin axis incli- 
nation and the orientation of the orbital plane. The 
planetary-mass companion 2M0122b has the only exo- 
planetary obliquity measurement to date and its obliq- 
uity may be high (Bryan et al. 2020b). A scenario that 
can produce this companion and give it a misaligned 
obliquity is formation via gravitational instability, where 
gravito-turbulence in the disk can torque fragment spin 
axes out of alignment (Bryan et al. 2020b). For the 
HR 8799 planets, formation via core accretion is pre- 
ferred based on occurrence rate statistics and atmo- 
spheric compositions (Konopacky et al. 2013; Barman 
et al. 2015; Nielsen et al. 2019; Molliére et al. 2020). 
In this case, obliquity excitation may arise from chaotic 
spin-orbit dynamics due to the near commensurability 
of nodal and spin precession frequencies (e.g. Neron de 
Surgy & Laskar 1997; Li & Batygin 2014; Shan & Li 
2018; Saillenfest et al. 2019). Specifically, consider the 
planet i = b,c,d,e, with mass mj, radius R;, spin fre- 
quency Q; = Q;/(Gm;/ R2)!/?, and semi-major axis aj, 
orbiting a host star of mass M,. Because the nodal pre- 
cession rate of planet i due to the other planets is of 
order (e.g. Pu & Lai 2018) 

" " 3Gm;a2. 

a 2 "p Vc 
where a< = min(a;,a;) and a> = max(a;, aj), while the 
spin-precession rate of planet i due to the host star's 
tidal torque is (assuming a fully convective planet, e.g. 
Lai 2014) 
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Figure 6. The final rotation rates relative to their breakup velocities as a function of mass for substellar objects. The left 
column assumes companion spins are aligned with their orbital planes, and the right column assumes randomly oriented spin 
axes. The top row only considers substellar companions excluding Jupiter and Saturn, the middle row includes Jupiter and 
Saturn, and the bottom row also includes free-floating substellar objects. In each plot, 50 possible power law relations using 
the fits from Table 4 are plotted as purple lines with their dispersions shaded in blue. HR 8799 d is green and e is yellow. 
HR 8799 c is in red, but we note that one should consider this an upper limit. Other substellar companions with robust spin 
measurements are plotted in black. Free-floating objects are plotted in gray. Rotational velocities derived from photometrically 
measured rotation periods are unchanged between the two plots. 


Table 3. Derived rotation rate values for the HR. 8799 planets. For non-detections of rotation, the 9596 upper limit or 596 lower limit is reported. 


Otherwise, the values reported follow the convention in Table 2. 
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Table 4. Power law fits to the rotation rate of gas giant planets. The values reported follow the convention in Table 
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and all spin frequencies have magnitudes 0; ~ 0.3, we 
have for the HR 8799 planets, wy/ay ~ 8.1, we/ae ~ 5.9, 
wa/da ~ 2.8, and we/Qe ~ 0.6, so all ratios of w;/a; are 
of order unity. Hence, this system may be susceptible to 
secular spin-orbit resonances (e.g. Millholland & Laugh- 
lin 2019; Bryan et al. 2020a). The coupled inclination 
and obliquity dynamics of the HR 8799 system will be 
investigated in a future work. 

In the following analysis, we considered two bounding 
assumptions: the spin axes being randomly oriented in 
space as was done in Bryan et al. (2020a), and spin axes 
being aligned with the orbital planes based on orbital 
inclination posteriors derived from the literature (Gin- 
ski et al. 2014; Bryan et al. 2016; Wang et al. 2018b; 
Pearce et al. 2019; Bryan et al. 2020b; Bowler et al. 
2020; Nowak, M. et al. 2020). We considered these 
to be bounding assumptions as reality likely does not 
have perfectly aligned obliquities nor a large fraction of 
retrograde obliquities. Two companions (SR 12 c and 
2M0249b) did not have measurements of orbital motion 
or photometrically derived rotation periods, so the pri- 
ors on their spin axes were assumed to be isotropic in 
both cases. Similarly, the free-floating objects were as- 
sumed to have isotropic spin axes orientations. 

In both cases, we computed the rotational velocities 
(v), with a preference of using photometrically derived 
rotation periods over vsin(i) measurements, since typ- 
ically radius uncertainties are smaller than the inclina- 
tion uncertainties. For the HR 8799 planets, we used an 
orbital inclination of i = 24? + 3° for the zero obliquity 
case (Wang et al. 2018b). We compared the rotation ve- 
locities against their break-up velocities, which we define 
simply as: 


Ubreak — GM/R, (12) 


where M is the mass of a object and R is its radius. This 
ignores the effect of oblateness, which will decrease the 
break-up velocity, with the exact decrease depending on 
the moment of inertia of the bodies (Marley & Sengupta 
2011; Tannock et al. 2021). To derive breakup velocities 
for the HR 8799 planets, we used masses of 7.20.7 Mjup 
for the inner three planets, a mass of 5.8+0.5 Mjup for 
HR 8799 b, and radii of 1.2 + 0.1 Ryup (Marois et al. 
2008, 2010; Wang et al. 2018b). The values of Upreak 
for each planet are listed in Table 3. Under both obliq- 
uity assumptions, we computed v, rotation period, and 
v/VUpreak, and listed the values for the HR 8799 plan- 
ets in Table 3. Note that we only reported 95% upper 
limits for HR. 8799 b and c. Although no rotation pe- 
riods have been measured from photometric monitoring 
of the HR 8799 planets (Apai et al. 2016; Biller et al. 
2021), our inferred rotation periods are > 3 hours for HR 
8799 d and e under both obliquity assumptions. This is 


consistent with the rotation periods of PSO J318.5-22 
and 2M1207b, the only objects similar in mass and age 
with photometrically measured rotation periods (Zhou 
et al. 2016; Biller et al. 2018). The rotation periods 
are also consistent with the picture that predicted mag- 
netic braking from the circumplanetary disk regulating 
the spin of the planet at very early times (< 1 Myr): 
Bryan et al. (2018, 2020a) found spins between 5-20% of 
breakup for planetary mass companions, Batygin (2018) 
predicted a terminal rotation period prediction of ~9 
hours, and Ginzburg & Chiang (2020) placed a maxi- 
mum rotation period of 29-43 hours. 

To compare these rotation rates against the whole 
population of substellar objects with measured rota- 
tion rates, we needed to account for the fact the ro- 
tation speed of an object depends on age. In our sam- 
ple, ages ranged from 108 to 10? years old. In the ab- 
sence of outside influence, an isolated body will spin 
up as its radius contracts due to radiative cooling with 
v c Ro} and v/vpreak X R-1/?, The current population 
of planetary mass objects has been shown to be consis- 
tent with this trend (Bryan et al. 2020a). We evolved 
the radii of all companions from their current radii to 
the radii predicted by hot-start evolutionary models at 
5 Gyr (Baraffe et al. 2003) and computed their final 
spin velocities relative to their final break-up velocities, 
Ufinal/Ubreak. The values of Vfnal/Vpreak for the HR. 8799 
planets are listed in Table 3. 

For the youngest objects in the sample, they may 
still host a circumplanetary disk (CPD) which combined 
with planetary magnetic fields is thought to regulate 
planetary spin (Batygin 2018; Ginzburg & Chiang 2020). 
CPD lifetimes are poorly constrained, so we do not know 
for sure which objects are currently experiencing mag- 
netic breaking. With evidence that CPDs and accretion 
occurs for the PDS 70 planets (Haffert et al. 2019; Isella 
et al. 2019), which are 8 + 1 Myr (Wang et al. 2020b), 
we assumed CPDs and regulation of planetary spin by 
magnetic braking occurs up to 10 Myr. For compan- 
ions younger than 10 Myr, we made one modification 
and assumed their rotational velocities are constant un- 
til 10 Myr. We note that magnetic braking is actually 
spinning down the planets, although the time depen- 
dence is weak (see Figure 1 of Ginzburg & Chiang 2020), 
so this is a rough approximation. 

We compared the vgnai/Upreak values of the HR 8799 
planets against the rest the gas giant and brown dwarf 
population in Figure 6 under our bounding obliquity as- 
sumptions and with three different subsets of the data 
(substellar companions excluding Jupiter and Saturn, 
substellar companions including Jupiter and Saturn, 
and all objects including free floating substellar objects). 
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Regardless of obliquity assumption and data subset, Fig- 
ure 6 appears to show a trend of increasing rotational 
velocity (relative to break-up velocity) with decreasing 
companion mass. A similar trend in v and mass for sub- 
stellar objects has been pointed out in previous works 
(Snellen et al. 2014; Scholz et al. 2018), but no significant 
trend was found when looking at 1-20 M),, planetary 
mass objects previously (Bryan et al. 2018; Zhou et al. 
2019; Xuan et al. 2020). Batygin (2018) and Ginzburg & 
Chiang (2020) suggested that lower mass planets may be 
less effective at ionizing their CPDs, making magnetic 
braking less effective in spinning down the planets and 
creating a spin speed that is mass dependent. To quan- 
tify the influence of mass on rotation rates, we followed 
the methodology from Wolfgang et al. (2016) where they 
used hierarchical Bayesian modeling to derive a proba- 
bilistic mass-radius relationship for transiting exoplan- 
ets. This approach accounts for uncertainties in both 
the independent and dependent variables by using their 
posteriors in the fit. This is important for our study as 
the constraints on mass and rotation speed can be weak 
and non-Gaussian, so point-estimates such as the me- 
dian and 6896 credible interval do not provide the full 
picture. We fit a model that is a power-law with intrinsic 
scatter of the following functional form: 


B 
Ufinal M 
Ubreak xd (r ( 10Mjup ) 7 k 
(13) 


Following the notation of Wolfgang et al. (2016), the ~ 
above indicates that the rotation rate is “drawn from 
the Normal distribution” with a mean (u) and standard 
deviation (c) that varies with companion mass (M). We 
fit for the power law constant (C), the power law index 
(3), and the fractional dispersion of the distribution at a 
given mass (øf). For both spin axis inclination assump- 
tions and for the three different data subsets, we esti- 
mated these population-level parameters and marginal- 
ized over the priors on rotation velocities and masses 
for the individual companions using the dynesty nested 
sampling package (Speagle 2020) with multiple bound- 
ing ellipsoids (Feroz et al. 2009) and random slice sam- 
pling (Neal 2003; Handley et al. 2015a,b). The priors 
and inferred values of the population-level parameters 
are listed in Table 4 and the fits excluding the Solar 
System gas giants are plotted in Figure 6. 

In all cases, we found a negative power law index was 
preferred. With the exception of the case where we 
looked at the substellar companions excluding Jupiter 
and Saturn and assuming zero obliquity, the 9596 cred- 
ible interval of 8 was entirely negative, although the 
mass dependence is weak: the 9596 credible interval for 


D is between -0.5 and 0 when including the Solar System 
gas giants, and between -1.2 and +0.1 when excluding 
them. This could be a tentative sign that the degree of 
ionization of the CPD by the planet itself affects the final 
rotation rate of a planet. However, there are other ways 
to create a mass dependence in v/Upreak: for example, 
the effect of oblateness would cause Vpreak to decrease by 
a factor that depends on the object’s moment of inertia, 
which itself is also mass dependent (Marley & Sengupta 
2011; Tannock et al. 2021). Even if it is true, the large 
scatter in the relation (c; ranging from 25-5096) indi- 
cates there would still be other variables at play. For 
example, since the planet is free to spin-up again after 
the CPD disperses in the magnetic braking scenario, dif- 
ferent CPD lifetimes could also contribute to the scatter 
in the population (Bryan et al. 2020a). Such a scatter 
is seen in young free-floating brown dwarfs where a ro- 
tation period dispersion was linked to the presence or 
absence of a disk at ~8 Myr ages (Scholz et al. 2018; 
Moore et al. 2019). In our fits, including the free-floating 
objects caused c; to reach the upper bound of the prior 
(5096 dispersion), which may indicate that mass is not a 
driving factor in shaping the spins of the free-floating 
objects, or that they are affected in a different way, 
since they formed through a different process than plan- 
etary companions. As the spins of more companions in 
the 1-10 Mjyp range are measured with instruments like 
KPIC, REACH, and HiRISE, we will be able to more 
confidently assess whether there is a trend in their ro- 
tation velocity with mass, and if there are other factors 
at play. 


6. CONCLUSIONS 


We obtained high spectral resolution (R ~ 35,000) K- 
band spectra of all four HR 8799 planets taken as part 
of the science-demonstration of KPIC, a new instrument 
that combines high contrast imaging with high resolu- 
tion spectroscopy to enable HDC techniques. By cross- 
correlating the spectra with molecular templates, we de- 
tected both CO and H20 at high significance (> 100 
combined) for HR 8799 c, d, and e, and tentatively de- 
tected HR 8799 b at 3c. These are the first detections 
of HR 8799 b, d, and e at high spectral resolution (R = 
10,000), where we can fully resolve the individual molec- 
ular lines in the planets’ atmospheres. The detection of 
HR 8799 c has an SNR that is 2x better while using 
3.5x less exposure time than the previous NIRSPAO 
detection (Wang et al. 2018a). With KPIC, we are able 
to access closer in planets such as HR 8799 e (385 mas 
from its star), which has previously never been detected 
at either medium or high resolution, demonstrating the 
HDC capabilities of KPIC. These are the most challeng- 
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ing directly imaged exoplanets characterized with high 
resolution spectroscopy to date given their flux ratios 
and angular separations: while 8 Pic b (Snellen et al. 
2014) is slightly closer in than HR, 8799 e, it is an order- 
of-magnitude brighter in the K band, which makes it 
easier to isolate from the diffracted stellar light. 

Beyond molecular detections via cross correlation, we 
developed a likelihood-based approach to jointly fit the 
stellar speckles and planet light in our high resolution 
data. We used the Bayes factor computed from our 
framework to assess detection probabilities that corrob- 
orated the molecular CCF analysis: definitive detections 
of HR 8799 c, d, and e, but a tentative detection of HR. 
8799 b where the No Planet model is still 9% as likely 
as the best model fit with a planet. Using our high 
resolution data alone, we inferred bulk properties from 
their atmospheres using the BT-Settl atmospheric mod- 
els (Allard et al. 2012). The measured radial velocities 
of the four planets are consistent with previous orbital 
constraints (Wang et al. 2018b; Ruffio et al. 2019). We 
found HR 8799 c, d, and e to have Teg constrained to 
between 1200 and 1600 K, which is somewhat consis- 
tent with previous works that used the B T-Settl mod- 
els (Bonnefoy et al. 2016; Greenbaum et al. 2018), but 
significantly higher than other literature fits that used 
different model grids. We also favored values of log(g) 
above 4.5 for HR. 8799 c and d, which are inconsistent 
with fits at low spectral resolution and predictions of 
mass and radius from evolutionary models. We specu- 
lated that deficiencies in the BT-Settl models could be 
partially responsible for our inferred parameters. Fu- 
ture work to study these planets with more accurate at- 
mospheric models would be able to determine whether 
our results are due to model or data systematics, but 
we were not able to perform such an analysis due to the 
limited number of publicly available forward model grids 
that both model clouds and support a spectral resolu- 
tion of 35,000. 

We found strong evidence for rotational broadening 
in our spectra of HR 8799 d and e, although more data 
is needed to make them decisive detections. Taken at 
face value, our inferred vsin(i) = 10.1*27 km/s for HR 
8799 d and vsin(i) = 15.0*22 km/s for HR 8799 e cor- 
respond to 0.2479 08 «c.i, and 0.36* 0 0s upreak assuming 
zero obliquities and [LIO hago and O17) a dbsenk 
assuming random oblitquities. We found a rotational 
velocity upper limit for HR 8799 c of < 14 km/s cor- 
responding to < 0.36vpreak in the zero obliquity case, 
but the spin of HR 8799 b was essentially unconstrained 
given the tentative detection. These spin measurements 
are consistent with the picture of magnetic braking reg- 
ulating the spins of gas giant planets during the planet 


formation process. When combining these rotation mea- 
surements with literature rotation measurements of sub- 
stellar companions as well as our Solar System gas gi- 
ants, we found a tentative trend of increasing rotation 
rate with decreasing planet mass that could point to 
magnetic braking being less efficient for lower mass plan- 
ets. 


6.1. Future Prospects 


More spin measurements of 1-10 Mjyp companions 
will test whether this mass-spin relation holds, or if there 
are other key parameters to control the spin rate of giant 
planets. In our analysis, we assumed two bounding as- 
sumptions on the obliquity of the planets to derive true 
rotation rates. Obliquity constraints on the HR 8799 
planets themselves may be possible in the future if their 
periods can be derived from rotational modulations in 
their light curves. This could be possible with new hard- 
ware designed to improve the photometric calibration of 
coronagraphic systems (Bos 2020). 

Future analysis work that moves beyond molecular de- 
tection and characterization of bulk parameters will pro- 
vide more insights into the HR 8799 planets. Retrievals 
of the chemical abundances of their atmospheres have 
been applied to broadband photometry and lower reso- 
lution spectra (Lavie et al. 2017; Molliére et al. 2020). 
The likelihood framework we constructed to fit the high 
spectral resolution KPIC data would be straightforward 
to include in these Bayesian retrievals, allowing us to 
make use of the information at all spectral resolutions 
in understanding the chemical makeup of these planets. 
KPIC is also being used to conduct a spectroscopic sur- 
vey of substellar companions which were previously too 
close to their bright host stars to study at high spectral 
resolution with slit spectrographs. Measuring planet 
characteristics like spin, orbital parameters, and com- 
position across this population will allow us to look for 
trends that would pinpoint the dominant processes in 
their formation and evolution. 

Lastly, future planned upgrades to KPIC (Pezzato 
et al. 2019; Jovanovic et al. 2020) are designed to im- 
prove the coupling of planet light, reduce the coupling 
of star light, and mitigate the instrument thermal back- 
ground. As we are limited by thermal background noise 
currently, increasing the amount of planet light reach- 
ing the detector and reducing the amount of thermal 
background photons will directly translate into improve- 
ments in the SNR of KPIC observations. Re-observing 
the HR 8799 planets with these upcoming upgrades 
would give us more precise planetary spin measurements 
and improved sensitivity to any trace methane abun- 
dances in their atmospheres to constrain the level of 
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disequilibrium chemistry. The initial phase of KPIC 
is the bare-bones hardware necessary to enable HDC 
techniques, and refined HDC techniques will allow for 
detailed characterization of all of the directly imaged 
exoplanets. 
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Figure 7. Corner plot for HR 8799 b of the posterior distribution of key parameters in the Full Model fit. 


APPENDIX 


A. CORNER PLOTS FOR FULL MODEL FITS 


We present corner plots for the Full Model fits for each planet: HR 8799 b in Figure 7, HR 8799 c in Figure 8, HR 
8799 d in Figure 9, and HR 8799 e in Figure 10. For each planet, we only showed six parameters, which are the main 
parameters of interest. These included the planet bulk properties Tig, log(g), planetary RV, vsin(i), and planetary 
flux. We also showed the measured speckle flux in Order 33. We did not show the other orders for simplicity as they 
have the same behavior. We also did not show the error inflation, offset, LSF size scaling, or linear coefficients c; for 
combining together the stellar spectra as they are nuisance parameters. The inclusion of all of these parameters would 
have limited the readability of these plots. 
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Figure 8. Corner plot for HR 8799 c of the posterior distribution of key parameters in the Full Model fit. 


B. FORWARD MODEL FIT AND RESIDUALS 


To better to assess the fits in detail, we created a version of Figure 4 that is zoomed into a 7 nm chunk of the 
Full Model fit in Order 33 (Figure 11). We can see good agreement between the data and forward model on a 
channel-to-channel basis. 

Another way to assess the fits are satisfactory is to look for any correlated structure in the residuals. This could be 
due to model mismatch or correlated noise sources that we had not accounted for. Although we did not visually see 
any structures in the residuals, the fact there are ~2000 data points per order means we can look for lower levels of 
correlated residual structure by computing the autocorrelation function (ACF): 

X; RiRi+e 
ACF(z) SURE (B1) 
Here, R is the 1D residual vector of length equal to the number of data points in the order, x is the lag in pixels (lag 
of zero results in ACF = 1), and 5°, iterates over all elements of R. For an infinitely long sequence of uncorrelated 
Gaussian noise, the ACF is the Kronecker delta function. Any correlated structure in the residuals would result in 
non-zero values in the wings of the ACF. 
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Figure 9. Corner plot for HR 8799 d of the posterior distribution of key parameters in the Full Model fit. 


From the fits in Section 4.3, we computed the fit residuals for the parameter set from the Full Model fit of each 
planet with highest likelihood. For each echelle order used, we computed the ACF of these residuals and plotted them 
in Figure 12. From visual inspection of the ACFs, we see that they are almost entirely consistent with a delta function. 
For HR 8799 e, we saw some power in the wings indicating small correlated structures in the residual (< 5% of the 
total scatter in the residuals). 
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10. Corner plot for HR. 8799 e of the posterior distribution of key parameters in the Full Model fit. 
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Figure 11. A version of Figure 4, except zoomed into a 7 nm chunk of the fit. The gaps in the data and models correspond to 
channels that were masked due to bad pixels. 
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Figure 12. ACF of the fit residuals for the best-fitting model for each planet. For each planet, all of the orders we fit to are 
plotted. The lack of discernible power in the wings of ACF for each order indicates we are dominated by uncorrelated noise. 


